Este documento contiene el código de R empleado para la elaboración del proyecto final de series de tiempo.
Los datos corresponden a la serie del Ãndice accionario Standard and Poor’s 500, y se obtuvieron del sÃmbolo GSPC de acuerdo a la información histórica disponible en https://finance.yahoo.com/quote/%5EGSPC. Se empleó el precio ajustado de cierre del Ãndice del 3 de enero de 2020 al 20 de abril de 2021.
Cargar las librerÃas requeridas:
library("readxl")
library("xts")
library("tseries")
library("forecast")
library("rugarch")
library("vars")
library("PerformanceAnalytics")
library("DMwR")
Importar datosy revisar datos:
datos <- read_excel("gspc.xlsx")
summary(datos)
## date close
## Min. :2000-01-03 00:00:00 Min. : 676.5
## 1st Qu.:2005-05-03 06:00:00 1st Qu.:1165.9
## Median :2010-08-26 12:00:00 Median :1394.8
## Mean :2010-08-27 12:04:34 Mean :1684.8
## 3rd Qu.:2015-12-21 18:00:00 3rd Qu.:2085.4
## Max. :2021-04-20 00:00:00 Max. :4185.5
Crear el objeto sp500 como una serie de tiempo (un objeto de la clase xts):
sp500 <- xts(x = datos$close, order.by = datos$date)
Verificar que la periodicidad de la serie sea correcta (se esperan datos diarios de 20 años):
periodicity(sp500)
## Daily periodicity from 2000-01-03 to 2021-04-20
Efectivamente, la serie de tiempo es de periodicidad diaria, con datos del 3 de enero de 2020 al 20 de abril de 2021.
Se grafica la serie de tiempo para una primera inspección visual:
plot(sp500)
Se prueba estacionariedad con la prueba aumentada de Dickey-Fuller:
adf.test(sp500, alternative = 's')
##
## Augmented Dickey-Fuller Test
##
## data: sp500
## Dickey-Fuller = -0.49676, Lag order = 17, p-value = 0.982
## alternative hypothesis: stationary
El valor \(p\) cercano a uno indica que no se puede rechazar la hipótesis nula de que la serie presenta raÃz unitaria.
Se calcula una serie que contiene las primeras diferencias para verificar si el proceso es integrado de orden uno. Esta serie se denomina sp500_1:
sp500_1 <- log(sp500) - log(lag(sp500))
sp500_1 <- na.omit(sp500_1)
Verificar periodicidad de la serie de primeras diferencias:
periodicity(sp500_1)
## Daily periodicity from 2000-01-04 to 2021-04-20
Se tiene una observación menos debido al cálculo de primeras diferencias.
Se repite la prueba aumentada de Dickey-Fuller con la serie de primeras diferencias:
adf.test(sp500_1, alternative = 's')
## Warning in adf.test(sp500_1, alternative = "s"): p-value smaller than printed p-
## value
##
## Augmented Dickey-Fuller Test
##
## data: sp500_1
## Dickey-Fuller = -17.641, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
El valor \(p\) menor a 0.01 indica que se rechaza la hipótesis nula de que la serie sp500_1 tiene raÃz unitaria. Como esta serie contiene las primeras diferencias de sp500, entonces esta última es una serie integrada de orden 1.
Gráfica de la serie de primeras diferencias:
plot(sp500_1)
Autocorrelación de la serie sp500:
sp500_acf <- acf(sp500)
Autocorrelación parcial de la serie sp500:
sp500_pacf <- pacf(sp500)
Autocorrelación de la serie sp500_1:
sp500_1_acf <- acf(sp500_1)
Autocorrelación parcial de la serie sp500_1:
sp500_1_pacf <- pacf(sp500_1)
Se realiza la prueba de Ljung-Box para deterninar si existe correlación serial en la serie estacionaria (sp500_1), para órdenes desde 1 hasta 30:
for (i in 1:30) {
prueba = Box.test(sp500_1, lag = i, type = "Ljung")
cat(i, "\t", round(prueba$statistic, 1), "\t", prueba$p.value, "\n")
}
## 1 68.8 1.110223e-16
## 2 68.9 1.110223e-15
## 3 70.1 3.996803e-15
## 4 73.3 4.551914e-15
## 5 74.3 1.276756e-14
## 6 81.9 1.44329e-15
## 7 87.8 3.330669e-16
## 8 92.5 1.110223e-16
## 9 100.3 0
## 10 100.3 0
## 11 101.5 1.110223e-16
## 12 113.8 0
## 13 115 0
## 14 115.1 0
## 15 135.7 0
## 16 166.5 0
## 17 166.7 0
## 18 176.3 0
## 19 176.3 0
## 20 176.3 0
## 21 178 0
## 22 178 0
## 23 178.2 0
## 24 178.4 0
## 25 178.4 0
## 26 179.6 0
## 27 184.9 0
## 28 185.1 0
## 29 185.6 0
## 30 186.6 0
Los valores \(p\) iguales a cero indican que se rechaza la hipótesis nula de que las observaciones no tienen correlación serial para órdenes menores o iguales a 30.
De acuerdo a la documentación del comando ARMA, la especificación del modelo ARMA(\(p\), \(q\)) en R es
\[ y_t = a_0 + a_1 y_{t-1} + \cdots + a_p y_{t-p} + b_1 e_{t-1} + \cdots + b_q e_{t-q} + e_t\]
Los comandos para los criterios de información en R sólo funcionan para estimaciones hechas mediante máxima verosimilitud (MLE), y la implementación del comando arma sólo es capaz de estimar usando mÃnimos cuadrados ordinarios (OLS). El comando que emplea MLE es Arima, por lo que es el que se utilizará para ajustar los modelos y obtener criterios de información. Dado que se empleará la serie sp500_1 (que tiene removida la raÃz unitaria presente en la serie original), el orden de integración del modelo debe ser cero (ARIMA(\(p\), \(0\), \(q\))), lo cual equivale a un modelo ARMA(\(p\), \(q\)).
Como primer ejercicio se estima el modelo ARMA(1, 1):
arima_1_0_1 <- Arima(sp500_1, order = c(1, 0, 1))
print(arima_1_0_1)
## Series: sp500_1
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.0071 -0.1225 2e-04
## s.e. 0.1140 0.1132 2e-04
##
## sigma^2 estimated as 0.0001546: log likelihood=15902.42
## AIC=-31796.83 AICc=-31796.82 BIC=-31770.49
Para determinar el mejor modelo empleando el criterio de información de Akaike (AIC), se estman 100 modelos ARMA(\(p\), \(q\)) diferentes, donde los parámetros de orden \(p\) y \(q\) varÃan independientemente en el rango de 1 hasta 10. Para cada modelo se obtiene el criterio de Akaike, y éstos se almacenan en la matriz criterios_aic, donde el número de fila corresponde al orden \(p\) y el número de columna corresponde al orden \(q\):
criterios_aic <- matrix(nrow = 10, ncol = 10)
for (p in 1:10) {
for (q in 1:10) {
modelo <- Arima(sp500_1, order = c(p, 0, q))
criterios_aic[p,q] <- modelo$aic
}
}
La matriz criterios_aic resulta ser:
print(criterios_aic)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -31796.83 -31794.85 -31797.92 -31799.09 -31797.62 -31805.47 -31804.77
## [2,] -31794.97 -31793.88 -31792.37 -31798.09 -31797.17 -31815.28 -31802.85
## [3,] -31793.83 -31791.83 -31793.82 -31801.75 -31799.75 -31815.01 -31811.51
## [4,] -31798.27 -31800.00 -31801.23 -31799.93 -31804.20 -31815.45 -31815.01
## [5,] -31797.20 -31797.85 -31799.79 -31805.14 -31824.08 -31822.41 -31819.64
## [6,] -31808.80 -31814.52 -31815.80 -31816.47 -31813.69 -31821.95 -31829.05
## [7,] -31807.15 -31804.92 -31811.22 -31815.32 -31815.08 -31829.58 -31826.89
## [8,] -31806.75 -31820.05 -31801.15 -31826.01 -31824.20 -31822.00 -31823.08
## [9,] -31808.56 -31830.77 -31825.19 -31820.85 -31824.14 -31833.60 -31832.09
## [10,] -31806.61 -31804.60 -31830.82 -31828.03 -31817.85 -31828.08 -31829.97
## [,8] [,9] [,10]
## [1,] -31803.20 -31804.44 -31802.50
## [2,] -31818.27 -31802.54 -31800.50
## [3,] -31815.39 -31817.91 -31817.64
## [4,] -31809.93 -31810.88 -31823.46
## [5,] -31821.80 -31830.96 -31841.80
## [6,] -31820.85 -31827.64 -31836.98
## [7,] -31835.50 -31830.42 -31834.89
## [8,] -31831.50 -31831.22 -31834.37
## [9,] -31831.05 -31844.89 -31840.76
## [10,] -31836.79 -31843.16 -31841.25
El valor mÃnimo es:
print(min(criterios_aic))
## [1] -31844.89
Que corresponde a la siguiente entrada:
which(criterios_aic == min(criterios_aic), arr.ind = TRUE)
## row col
## [1,] 9 9
Por lo tanto, el modelo óptimo es el ARMA(9, 9). Se ajusta dicho modelo a los datos:
arima_9_0_9 <- Arima(sp500_1, order = c(9, 0, 9))
print(arima_9_0_9)
## Series: sp500_1
## ARIMA(9,0,9) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## 0.1031 -0.3081 0.5003 0.3950 -0.2157 -0.5844 -0.1778 -0.4054
## s.e. 0.0706 0.0387 0.0796 0.0872 0.0602 0.0674 0.0805 0.0611
## ar9 ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8
## 0.3980 -0.2165 0.3167 -0.5186 -0.3737 0.2616 0.5165 0.1524 0.3641
## s.e. 0.0694 0.0690 0.0623 0.0937 0.1013 0.0786 0.0806 0.0754 0.0649
## ma9 mean
## -0.3697 2e-04
## s.e. 0.0747 1e-04
##
## sigma^2 estimated as 0.0001528: log likelihood=15942.44
## AIC=-31844.89 AICc=-31844.73 BIC=-31713.17
residuales <- checkresiduals(arima_9_0_9)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(9,0,9) with non-zero mean
## Q* = 21.2, df = 3, p-value = 9.567e-05
##
## Model df: 19. Total lags used: 22
Valor \(p\) de la prueba de Ljung-Box sobre los residuales del modelo ajustado:
print(residuales$p.value)
## [1] 9.567302e-05
Debido a las que presenta la serie con periodicidad diaria, se convierte a datos mensuales empleando el valor de cierre (del último dÃa del mes) como valor mensual. Esta nueva serie se denomina sp500m:
sp500m <- to.monthly(sp500)
sp500m <- sp500m$sp500.Close
Gráfica de la serie mensual:
plot(sp500m)
Prueba de Dickey-Fuller aumentada:
adf.test(sp500m, alternative = 's')
## Warning in adf.test(sp500m, alternative = "s"): p-value greater than printed p-
## value
##
## Augmented Dickey-Fuller Test
##
## data: sp500m
## Dickey-Fuller = 0.24955, Lag order = 6, p-value = 0.99
## alternative hypothesis: stationary
No se puede rechazar la hipótesis nula de que la serie tiene raÃz unitaria. Se obtiene la serie sp500m_1 que contiene las primeras diferencias:
sp500m_1 <- log(sp500m) - log(lag(sp500m))
sp500m_1 <- na.omit(sp500m_1)
Gráfica de sp500m_1:
plot(sp500m_1)
Prueba de DIckey-Fuller aumentada para sp500m_1:
adf.test(sp500m_1)
## Warning in adf.test(sp500m_1): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: sp500m_1
## Dickey-Fuller = -5.6307, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
Se rechaza la hipótesis nula de que la serie sp500m_1 tiene raÃz unitaria.
Función de autocorrelación de sp500m_1:
sp500m_1_acf <- acf(sp500m_1)
Función de autocorrelación parcial de sp500m_1:
sp500m_1_pacf <- pacf(sp500m_1)
Prueba de Ljung-Box para determinar si existe autocorrelación serial:
for (i in 1:30) {
prueba = Box.test(sp500m_1, lag = i, type = "Ljung")
cat(i, "\t", round(prueba$statistic, 1), "\t", prueba$p.value, "\n")
}
## 1 1.5 0.2157794
## 2 2.8 0.2453356
## 3 4.4 0.2173632
## 4 5.7 0.223933
## 5 5.9 0.3124418
## 6 7.5 0.2761677
## 7 9.1 0.2486254
## 8 9.6 0.2921704
## 9 11.5 0.2440293
## 10 11.6 0.3143923
## 11 11.6 0.3937488
## 12 11.7 0.4724017
## 13 12.5 0.4898004
## 14 13.4 0.4972188
## 15 15.8 0.3935624
## 16 16.9 0.3952379
## 17 17 0.4555091
## 18 17.6 0.484824
## 19 21.3 0.319761
## 20 23 0.287934
## 21 23.2 0.3346201
## 22 23.7 0.3619937
## 23 23.7 0.4187952
## 24 24.1 0.4574857
## 25 26.3 0.3933144
## 26 28.8 0.3217141
## 27 29.2 0.3516093
## 28 29.8 0.3710967
## 29 30.4 0.392367
## 30 30.4 0.4432944
Se repite la búsqueda del modelo ARMA óptimo para la nueva serie sp500m_1, esta vez variando los parámetros \(p\) y \(q\) entre 1 y 10:
criterios_aic_m <- matrix(nrow = 10, ncol = 10)
for (p in 1:10) {
for (q in 1:10) {
modelo <- Arima(sp500m_1, order = c(p, 0, q), method = 'ML')
criterios_aic_m[p,q] <- modelo$aic
}
}
La matriz criterios_aic_m resulta ser:
print(criterios_aic_m)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -867.2541 -865.3633 -864.7100 -862.7104 -861.8807 -860.7557 -859.5645
## [2,] -865.3693 -863.3677 -862.7102 -864.6792 -866.0731 -864.1559 -861.7982
## [3,] -864.6408 -866.7266 -860.7138 -862.9615 -864.1541 -862.1545 -860.1897
## [4,] -862.7110 -864.7479 -863.3447 -866.0930 -864.2054 -860.1826 -863.0429
## [5,] -862.2362 -865.7666 -865.8954 -864.2013 -862.1549 -858.2120 -857.1438
## [6,] -861.5606 -863.7697 -861.8072 -859.7982 -858.1328 -856.1995 -860.2383
## [7,] -860.7160 -861.5805 -859.8537 -860.7448 -860.7414 -859.1952 -862.6002
## [8,] -859.7129 -860.6694 -858.0276 -858.8732 -858.1847 -862.6773 -855.7247
## [9,] -858.6157 -856.6158 -854.6661 -859.8558 -857.7075 -860.8289 -858.8275
## [10,] -856.6165 -854.6157 -856.7950 -855.3184 -854.8144 -858.8293 -856.2785
## [,8] [,9] [,10]
## [1,] -859.6983 -859.3776 -857.9054
## [2,] -860.6344 -858.6660 -864.0275
## [3,] -858.2261 -856.6941 -861.3563
## [4,] -860.0398 -854.6932 -854.4539
## [5,] -858.3866 -857.8570 -855.7522
## [6,] -860.4172 -860.8915 -858.9088
## [7,] -857.2897 -858.8627 -856.9329
## [8,] -857.2197 -857.0902 -855.2309
## [9,] -857.1309 -857.0570 -856.0372
## [10,] -854.7100 -855.8460 -859.5798
El valor mÃnimo es:
print(min(criterios_aic_m))
## [1] -867.2541
Que corresponde a la siguiente entrada:
which(criterios_aic_m == min(criterios_aic_m), arr.ind = TRUE)
## row col
## [1,] 1 1
Por lo tanto, el modelo óptimo es el ARMA(1, 1). Se ajusta dicho modelo a los datos:
arima_m_1_0_1 <- Arima(sp500m_1, order = c(1, 0, 1))
print(arima_m_1_0_1)
## Series: sp500m_1
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## -0.6531 0.7552 0.0043
## s.e. 0.2674 0.2323 0.0029
##
## sigma^2 estimated as 0.001914: log likelihood=437.63
## AIC=-867.25 AICc=-867.09 BIC=-853.09
Residuales:
residuales_m <- checkresiduals(arima_m_1_0_1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 22.279, df = 21, p-value = 0.3836
##
## Model df: 3. Total lags used: 24
Valor \(p\) de la prueba de Ljung-Box sobre los residuales del modelo ajustado:
print(residuales_m$p.value)
## [1] 0.3835919
Por tanto los residuales parecen ser ruido blanco.
Pronóstico dentro de la muestra:
training_set <- sp500m_1['/2018-12-31']
evaluation_set <- sp500m_1['2019-01-01/']
training_model <- Arima(training_set, order = c(1, 0, 1))
pronostico <- xts(predict(training_model, length(evaluation_set))$pred, order.by = index(evaluation_set))
pronostico_in <- merge(pronostico, evaluation_set, join='inner')
plot.xts(pronostico_in, screens = factor(1, 1), legend.loc = 'bottomright')
Medidas de error del pronóstico dentro de la muestra:
medidas_error = regr.eval(coredata(evaluation_set), coredata(predict(training_model, length(evaluation_set))$pred))
print(round(medidas_error, 4))
## mae mse rmse mape
## 0.0449 0.0031 0.0561 1.0233
Pronóstico fuera de la muestra:
pronostico <- xts(predict(training_model, 12)$pred, order.by = seq(as.Date("2021-06-01"),length=12,by="months"))
pronostico_out <- merge(pronostico, sp500_1, join = 'outer', fill = 0)['2021-01-01/']
plot.xts(pronostico_out, screens = factor(1, 1), legend.loc = 'bottomright')
Se estiman 25 modelos GARCH diferentes variando los parámetros \(p\) y \(q\) entre 1 y 5. Para cada modelo se calcula el criterio de información bayesiano (BIC) y éstos se almacenan en la matriz criterios_bic, donde el número de fila corresponde al orden de \(p\) y el número de columna corresponde al orden de \(q\):
criterios_bic = matrix(nrow = 5, ncol= 5)
for (p in 1:5) {
for (q in 1:5) {
# Especificación del modelo:
garch_spec <- ugarchspec(
variance.model = list(
model = "sGARCH",
garchOrder = c(p, q)),
mean.model = list(
armaOrder = c(p, q)),
distribution.model = "norm")
# Ajuste del modelo:
garch_model <- ugarchfit(spec = garch_spec, data = sp500_1)
# Cálculo del BIC:
criterios_bic[p, q] = infocriteria(garch_model)[2]
}
}
Matriz de criterios BIC:
criterios_bic
## [,1] [,2] [,3] [,4] [,5]
## [1,] -6.449218 -6.446946 -6.443827 -6.440720 -6.437653
## [2,] -6.449317 -6.446842 -6.443833 -6.440217 -6.437456
## [3,] -6.445972 -6.443891 -6.440953 -6.439616 -6.434335
## [4,] -6.442883 -6.440892 -6.437719 -6.435750 -6.432870
## [5,] -6.439802 -6.437584 -6.438893 -6.432912 -6.430124
El valor mÃnimo es:
print(min(criterios_bic))
## [1] -6.449317
Que corresponde a la siguiente entrada:
which(criterios_bic == min(criterios_bic), arr.ind = TRUE)
## row col
## [1,] 2 1
Por lo que el modelo óptimo es el GARCH(2, 1). Este modelo se ajusta en el objeto garch_2_1:
garch_spec <- ugarchspec(
variance.model = list(
model = "sGARCH",
garchOrder = c(2, 1)),
mean.model = list(
armaOrder = c(2, 1)),
distribution.model = "norm")
garch_2_1 <- ugarchfit(spec = garch_spec, data = sp500_1)
Resultado del ajuste:
show(garch_2_1)
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(2,1)
## Mean Model : ARFIMA(2,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000621 0.000075 8.3143 0.000000
## ar1 0.891415 0.014314 62.2750 0.000000
## ar2 0.043317 0.014708 2.9452 0.003228
## ma1 -0.956180 0.001249 -765.3050 0.000000
## omega 0.000003 0.000001 3.3912 0.000696
## alpha1 0.081623 0.013697 5.9591 0.000000
## alpha2 0.064236 0.016971 3.7849 0.000154
## beta1 0.834458 0.012874 64.8185 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000621 0.000076 8.1919 0.000000
## ar1 0.891415 0.013650 65.3032 0.000000
## ar2 0.043317 0.014219 3.0464 0.002316
## ma1 -0.956180 0.000397 -2408.4911 0.000000
## omega 0.000003 0.000004 0.6531 0.513690
## alpha1 0.081623 0.024367 3.3498 0.000809
## alpha2 0.064236 0.035469 1.8111 0.070133
## beta1 0.834458 0.044652 18.6880 0.000000
##
## LogLikelihood : 17308.84
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.4592
## Bayes -6.4493
## Shibata -6.4592
## Hannan-Quinn -6.4557
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.072 0.3005
## Lag[2*(p+q)+(p+q)-1][8] 4.051 0.7620
## Lag[4*(p+q)+(p+q)-1][14] 7.657 0.4336
## d.o.f=3
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.211 0.2712
## Lag[2*(p+q)+(p+q)-1][8] 3.595 0.5770
## Lag[4*(p+q)+(p+q)-1][14] 7.396 0.4557
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.5986 0.500 2.000 0.4391
## ARCH Lag[6] 2.0371 1.461 1.711 0.4826
## ARCH Lag[8] 2.8568 2.368 1.583 0.5706
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 12.6708
## Individual Statistics:
## mu 0.70986
## ar1 0.12406
## ar2 0.09757
## ma1 0.18289
## omega 0.71000
## alpha1 0.19492
## alpha2 0.56190
## beta1 0.78566
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.89 2.11 2.59
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 3.2301 1.245e-03 ***
## Negative Sign Bias 0.2194 8.264e-01
## Positive Sign Bias 2.3019 2.138e-02 **
## Joint Effect 39.5560 1.323e-08 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 196.6 1.638e-31
## 2 30 226.0 2.143e-32
## 3 40 259.8 2.035e-34
## 4 50 298.9 1.452e-37
##
##
## Elapsed time : 0.95836
Se observa que todos los coeficientes son significativos.
Se analiza qué sucede al aumentar el número de parámetros a \(p=3\) y \(q=2\):
garch_spec <- ugarchspec(
variance.model = list(
model = "sGARCH",
garchOrder = c(3, 2)),
mean.model = list(
armaOrder = c(3, 2)),
distribution.model = "norm")
garch_3_2 <- ugarchfit(spec = garch_spec, data = sp500_1)
show(garch_3_2)
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(3,2)
## Mean Model : ARFIMA(3,0,2)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000621 0.000072 8.56391 0.000000
## ar1 0.594759 0.014221 41.82339 0.000000
## ar2 0.304970 0.016806 18.14609 0.000000
## ar3 0.016698 0.014684 1.13714 0.255482
## ma1 -0.659325 0.001306 -504.96046 0.000000
## ma2 -0.285140 0.003228 -88.34446 0.000000
## omega 0.000005 0.000001 3.99799 0.000064
## alpha1 0.079003 0.013758 5.74248 0.000000
## alpha2 0.143467 0.013550 10.58773 0.000000
## alpha3 0.025565 0.017609 1.45182 0.146552
## beta1 0.069598 0.092797 0.75001 0.453251
## beta2 0.649545 0.087091 7.45821 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000621 0.000075 8.29435 0.000000
## ar1 0.594759 0.013379 44.45413 0.000000
## ar2 0.304970 0.017513 17.41380 0.000000
## ar3 0.016698 0.014742 1.13269 0.257345
## ma1 -0.659325 0.001295 -509.27693 0.000000
## ma2 -0.285140 0.000836 -341.14302 0.000000
## omega 0.000005 0.000005 1.00796 0.313473
## alpha1 0.079003 0.025352 3.11620 0.001832
## alpha2 0.143467 0.018735 7.65778 0.000000
## alpha3 0.025565 0.020267 1.26141 0.207160
## beta1 0.069598 0.138611 0.50211 0.615590
## beta2 0.649545 0.164953 3.93775 0.000082
##
## LogLikelihood : 17311.48
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.4586
## Bayes -6.4439
## Shibata -6.4587
## Hannan-Quinn -6.4535
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.084 0.2977
## Lag[2*(p+q)+(p+q)-1][14] 7.572 0.4419
## Lag[4*(p+q)+(p+q)-1][24] 15.389 0.1388
## d.o.f=5
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.638 0.2006
## Lag[2*(p+q)+(p+q)-1][14] 5.989 0.6404
## Lag[4*(p+q)+(p+q)-1][24] 9.424 0.7595
## d.o.f=5
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[6] 1.356 0.500 2.000 0.2443
## ARCH Lag[8] 1.851 1.480 1.774 0.5471
## ARCH Lag[10] 5.696 2.424 1.650 0.2136
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 6.2489
## Individual Statistics:
## mu 0.66619
## ar1 0.13507
## ar2 0.06704
## ar3 0.21539
## ma1 0.19551
## ma2 0.13707
## omega 0.17338
## alpha1 0.14043
## alpha2 0.42789
## alpha3 0.74637
## beta1 0.77642
## beta2 0.82597
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 2.69 2.96 3.51
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 3.2132 1.320e-03 ***
## Negative Sign Bias 0.3199 7.491e-01
## Positive Sign Bias 2.2909 2.201e-02 **
## Joint Effect 39.9379 1.098e-08 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 197.3 1.184e-31
## 2 30 230.7 2.673e-33
## 3 40 258.5 3.601e-34
## 4 50 286.6 2.598e-35
##
##
## Elapsed time : 1.314727
Normalidad de los residuales:
plot(garch_2_1, which = 9)
Los residuales no parecen seguir una distribución normal.
Autocorrelación de los residuales:
plot(garch_2_1, which = 10)
Gráfico de la volatilidad con 2 desviaciones estándar:
plot(garch_2_1, which = 1)
El valor en riesgo de la serie es (para \(p\) igual a 0.90, 0.95 y 0.99):
for (prob in c(0.90, 0.95, 0.99)) {
print(c(prob, VaR(sp500_1, p = prob, method = 'historical')))
}
## [1] 0.90000000 -0.01301216
## [1] 0.95000000 -0.01913197
## [1] 0.99000000 -0.03520641
Serie del retorno de largo plazo del Tesoro de EUA (compuesto a 10 años): https://www.treasury.gov/resource-center/data-chart-center/interest-rates/pages/TextView.aspx?data=longtermrateAll
Importar datos y revisar datos:
datos_ltc <- read_excel("lt_composite.xlsx")
summary(datos_ltc)
## date lt_composite
## Min. :2000-01-03 00:00:00 Min. :0.970
## 1st Qu.:2005-04-28 18:00:00 1st Qu.:2.740
## Median :2010-08-16 12:00:00 Median :3.810
## Mean :2010-08-17 03:47:58 Mean :3.798
## 3rd Qu.:2015-12-07 06:00:00 3rd Qu.:4.820
## Max. :2021-03-31 00:00:00 Max. :6.910
Crear el objeto ltc como una serie de tiempo (un objeto de la clase xts):
ltc <- xts(x = datos_ltc$lt_composite, order.by = datos_ltc$date)
Verificar que la periodicidad de la serie sea correcta (se esperan datos diarios de 20 años):
periodicity(ltc)
## Daily periodicity from 2000-01-03 to 2021-03-31
Gráfica de la serie de tiempo:
plot(ltc)
Se prueba estacionariedad con la prueba ADF:
adf.test(ltc, alternative = 's')
## Warning in adf.test(ltc, alternative = "s"): p-value smaller than printed p-
## value
##
## Augmented Dickey-Fuller Test
##
## data: ltc
## Dickey-Fuller = -4.3485, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
El valor \(p\) de 0.01 indica que se rechaza la hipótesis nula de que la serie presenta raÃz unitaria. Aun asÃ, para que esta serie sea consistente con la sp500_1, se laculan las primeras diferencias en una nueva serie ltc_1:
ltc_1 <- ltc - lag(ltc)
ltc_1 <- na.omit(ltc_1)
Verificar periodicidad de la serie de primeras diferencias:
periodicity(ltc_1)
## Daily periodicity from 2000-01-04 to 2021-03-31
Se tiene una observación menos debido al cálculo de primeras diferencias.
Se repite la prueba aumentada de Dickey-Fuller con la serie de primeras diferencias:
adf.test(ltc_1)
## Warning in adf.test(ltc_1): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ltc_1
## Dickey-Fuller = -15.967, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
La serie conserva su estacionariedad original.
Gráfico de la serie de primeras diferencias:
plot(ltc_1)
Se encuentra el inicio y el final comunes de ambas series:
inicio <- max(min(index(sp500_1)), min(index(ltc_1)))
final <- min(max(index(sp500_1)), max(index(ltc_1)))
Estas fechas son, respectivamente:
print(c(inicio, final))
## [1] "2000-01-03 18:00:00 CST" "2021-03-30 18:00:00 CST"
res <- merge(sp500, ltc, join='inner')
VARselect(res)
## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 10 10 10 10
##
## $criteria
## 1 2 3 4 5 6 7
## AIC(n) 0.1771518 0.1528343 0.1477548 0.1478244 0.1460910 0.1473937 0.1343622
## HQ(n) 0.1797549 0.1571728 0.1538286 0.1556336 0.1556355 0.1586736 0.1473775
## SC(n) 0.1846005 0.1652488 0.1651351 0.1701705 0.1734028 0.1796713 0.1716056
## FPE(n) 1.1938123 1.1651320 1.1592287 1.1593094 1.1573015 1.1588101 1.1438071
## 8 9 10
## AIC(n) 0.1229181 0.1187283 0.1129452
## HQ(n) 0.1376687 0.1352143 0.1311666
## SC(n) 0.1651273 0.1659033 0.1650860
## FPE(n) 1.1307918 1.1260640 1.1195707
var_10 <- VAR(res, p = 10, type = 'const')
summary(var_10)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: sp500, ltc
## Deterministic variables: const
## Sample size: 5296
## Log Likelihood: -15286.476
## Roots of the characteristic polynomial:
## 1 0.9975 0.8676 0.8676 0.7582 0.7582 0.755 0.755 0.7129 0.6897 0.6897 0.6869 0.6869 0.6728 0.6728 0.6495 0.6495 0.6369 0.6369 0.01075
## Call:
## VAR(y = res, p = 10, type = "const")
##
##
## Estimation results for equation sp500:
## ======================================
## sp500 = sp500.l1 + ltc.l1 + sp500.l2 + ltc.l2 + sp500.l3 + ltc.l3 + sp500.l4 + ltc.l4 + sp500.l5 + ltc.l5 + sp500.l6 + ltc.l6 + sp500.l7 + ltc.l7 + sp500.l8 + ltc.l8 + sp500.l9 + ltc.l9 + sp500.l10 + ltc.l10 + const
##
## Estimate Std. Error t value Pr(>|t|)
## sp500.l1 0.87675 0.01448 60.543 < 2e-16 ***
## ltc.l1 -1.87876 5.44382 -0.345 0.730020
## sp500.l2 0.17905 0.01933 9.261 < 2e-16 ***
## ltc.l2 -9.25540 7.72063 -1.199 0.230664
## sp500.l3 -0.02093 0.01935 -1.081 0.279546
## ltc.l3 7.64278 7.71843 0.990 0.322122
## sp500.l4 -0.08706 0.01915 -4.545 5.61e-06 ***
## ltc.l4 3.69813 7.71172 0.480 0.631569
## sp500.l5 0.04667 0.01914 2.439 0.014778 *
## ltc.l5 -4.34354 7.70504 -0.564 0.572964
## sp500.l6 -0.09209 0.01915 -4.810 1.55e-06 ***
## ltc.l6 16.75382 7.70438 2.175 0.029706 *
## sp500.l7 0.19857 0.01916 10.366 < 2e-16 ***
## ltc.l7 -28.94440 7.70762 -3.755 0.000175 ***
## sp500.l8 -0.15877 0.01939 -8.190 3.24e-16 ***
## ltc.l8 18.66746 7.71211 2.421 0.015531 *
## sp500.l9 0.14443 0.01936 7.460 1.01e-13 ***
## ltc.l9 -6.24425 7.71494 -0.809 0.418338
## sp500.l10 -0.08713 0.01450 -6.011 1.97e-09 ***
## ltc.l10 2.99994 5.43673 0.552 0.581115
## const 4.75519 1.98274 2.398 0.016506 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 20.41 on 5275 degrees of freedom
## Multiple R-Squared: 0.9992, Adjusted R-squared: 0.9992
## F-statistic: 3.208e+05 on 20 and 5275 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation ltc:
## ====================================
## ltc = sp500.l1 + ltc.l1 + sp500.l2 + ltc.l2 + sp500.l3 + ltc.l3 + sp500.l4 + ltc.l4 + sp500.l5 + ltc.l5 + sp500.l6 + ltc.l6 + sp500.l7 + ltc.l7 + sp500.l8 + ltc.l8 + sp500.l9 + ltc.l9 + sp500.l10 + ltc.l10 + const
##
## Estimate Std. Error t value Pr(>|t|)
## sp500.l1 -1.181e-04 3.864e-05 -3.057 0.00225 **
## ltc.l1 1.005e+00 1.453e-02 69.183 < 2e-16 ***
## sp500.l2 1.292e-04 5.159e-05 2.504 0.01230 *
## ltc.l2 -5.484e-02 2.060e-02 -2.662 0.00779 **
## sp500.l3 -4.919e-05 5.164e-05 -0.952 0.34089
## ltc.l3 4.448e-02 2.060e-02 2.160 0.03085 *
## sp500.l4 -4.633e-05 5.111e-05 -0.906 0.36471
## ltc.l4 7.148e-04 2.058e-02 0.035 0.97229
## sp500.l5 8.841e-05 5.107e-05 1.731 0.08347 .
## ltc.l5 -1.864e-03 2.056e-02 -0.091 0.92778
## sp500.l6 3.611e-05 5.109e-05 0.707 0.47967
## ltc.l6 5.426e-03 2.056e-02 0.264 0.79186
## sp500.l7 -2.949e-05 5.112e-05 -0.577 0.56408
## ltc.l7 1.737e-03 2.057e-02 0.084 0.93271
## sp500.l8 -1.946e-05 5.173e-05 -0.376 0.70680
## ltc.l8 -3.623e-02 2.058e-02 -1.761 0.07835 .
## sp500.l9 1.025e-04 5.166e-05 1.984 0.04730 *
## ltc.l9 3.062e-02 2.059e-02 1.487 0.13703
## sp500.l10 -9.559e-05 3.868e-05 -2.471 0.01349 *
## ltc.l10 2.915e-03 1.451e-02 0.201 0.84076
## const 1.021e-02 5.291e-03 1.931 0.05359 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.05447 on 5275 degrees of freedom
## Multiple R-Squared: 0.9982, Adjusted R-squared: 0.9982
## F-statistic: 1.462e+05 on 20 and 5275 DF, p-value: < 2.2e-16
##
##
##
## Covariance matrix of residuals:
## sp500 ltc
## sp500 416.6996 0.354484
## ltc 0.3545 0.002967
##
## Correlation matrix of residuals:
## sp500 ltc
## sp500 1.0000 0.3188
## ltc 0.3188 1.0000
plot(var_10)
arch(var_10)
## Warning: Function 'arch' is deprecated; use 'arch.test' instead.
## See help("vars-deprecated") and help("arch-deprecated") for more information.
##
## ARCH (multivariate)
##
## data: Residuals of VAR object x
## Chi-squared = 3086.2, df = 45, p-value < 2.2e-16
ir <- irf(var_10, seed = 1)
plot(ir)
ir$runs
## [1] 100
ir$irf
## $sp500
## sp500 ltc
## [1,] 20.41322 0.01736539
## [2,] 17.86466 0.01504078
## [3,] 19.12880 0.01469069
## [4,] 19.50848 0.01375604
## [5,] 18.39539 0.01204305
## [6,] 18.56075 0.01236441
## [7,] 16.69747 0.01310932
## [8,] 18.86988 0.01347808
## [9,] 17.06684 0.01249074
## [10,] 19.27205 0.01472273
## [11,] 18.90516 0.01422970
##
## $ltc
## sp500 ltc
## [1,] 0.00000000 0.05162921
## [2,] -0.09699875 0.05188628
## [3,] -0.66037459 0.04932462
## [4,] -0.77465724 0.04908652
## [5,] -0.75664287 0.04898153
## [6,] -0.98155554 0.04869442
## [7,] -0.27452383 0.04873090
## [8,] -1.18072949 0.04871283
## [9,] -0.90665051 0.04697710
## [10,] -1.15325776 0.04664313
## [11,] -1.19386703 0.04671415